#notes for MK tests

chick.seqs <- read.BStringSet("/Paterson/Datafiles/grouse/seq_cap/blast_out/chick_transcripts.fasta")
chick.prot.cDNA <- read.table("/Paterson/Datafiles/grouse/seq_cap/blast_out/chick_prot_cDNA.txt",header=T,stringsAsFactors=F)




getCDS <- function(dir.embl="/Paterson/Datafiles/grouse/embl_exon2",seq.name,single.cds=FALSE,file.name=NULL){
	#AllSNPref is an object (dataframe) from baseCoverage
	#expects embl files in dir.embl, perhaps created by add.CDS
		
	
	if(is.null(file.name)) {
		#fasta.file.name<-paste(dir.fasta,"/",seq.name,"_matched.fasta",sep="")
		embl.file.name<-paste(dir.embl,"/",seq.name,"_exon.embl",sep="")
		}else{
		#fasta.file.name <- file.name[1]
		embl.file.name <- file.name[1]
		}
	
	if(!file.exists(embl.file.name)){
		warning(embl.file.name, " does not exist... returning\n")
		return(c(NA,NA))
		}

	#require(Biostrings)
	
	#for some reason won't read straight into buffer
	
	#ref.seq <- readEmblSeq(embl.file.name)
	#ref.seq.set <- read.DNAStringSet(fasta.file.name)
	#ref.seq <- DNAString(ref.seq)
	
	#extract feature table
	tmp.features <- system(paste("grep \'FT\'",embl.file.name,"|sed \'s/FT[ \t]//\'"),intern=TRUE)
	cds.lines <- grep('^[ ]*CDS ',tmp.features)
	product.lines <- grep('/product',tmp.features)
	grep.cds <- character(length=length(cds.lines))
	for(cdsi in 1:length(grep.cds)){
		grep.cds[cdsi] <- paste(tmp.features[cds.lines[cdsi]:(product.lines[cdsi]-1)],collapse="")
		
		}
	#extract coding seqs
	grep.cds <- sub('^[ ]*CDS ','',grep.cds)
	grep.cds <- gsub('[ ]*','',grep.cds)

	#grep.cds <- system(paste("grep \'CDS\'",embl.file.name,"|sed \'s/FT[ \t]*CDS[ \t]*//\'"),intern=TRUE)
	complementCDS <- FALSE
	if(substr(grep.cds[1],1,4)=="comp"){
		complementCDS <- TRUE
		}

	ref.seq <- readEmblSeq(embl.file.name)
	
	grep.cds <- sub("complement\\(","",grep.cds)
	grep.cds <- sub("join\\(","",grep.cds)
	grep.cds <- gsub(")","",grep.cds)

	#find masks
	mask.lines <- grep("/note=\"mask\"",tmp.features) -1
	if(length(mask.lines)>0){
		grep.mask <- tmp.features[mask.lines]
		grep.mask <- sub('^[ ]*misc_feature','',grep.mask)
		grep.mask <- gsub('[ ]*','',grep.mask)
		
		mask.list <- strsplit(grep.mask,"\\.\\.")
		mask.list <- lapply(mask.list,as.numeric)
		mask.list <- lapply(X=mask.list,FUN=function(X){
			if(length(X)==1){
				return(c(X,X))
				}else{
				return(X)
				}
			})
		mask.list <- matrix(unlist(mask.list),ncol=2,byrow=T)
		colnames(mask.list) <- c("start","end")
		mask.list <- as.data.frame(mask.list) 
		#output starts and ends of mask as a dataframe
		tmp.loc <- numeric(0)
		for(i in 1:nrow(mask.list)) tmp.loc <- c(tmp.loc,seq(mask.list[i,1],mask.list[i,2]))
		for(i in 1:length(tmp.loc)){
			ref.seq <- paste(substr(ref.seq,1,(tmp.loc[i]-1)),"#",substr(ref.seq,(tmp.loc[i]+1),nchar(ref.seq)),sep="")
			#will break if masked nucl at end of sequence
			}

		}

	#from mutateSeq
	
	
	tmp.from <- list()
	tmp.to <- list()
	for(i in 1:length(grep.cds)){
		tmp.from[[i]] <- as.numeric(sapply(X=strsplit(grep.cds[i],","),FUN=function(X){sub('\\.\\.[0-9]*$','',X)}))
		tmp.from[[i]] <- sort(tmp.from[[i]])
		tmp.to[[i]] <- as.numeric(sapply(X=strsplit(grep.cds[i],","),FUN=function(X){sub('^[0-9]*\\.\\.','',X)}))	
		tmp.to[[i]] <- sort(tmp.to[[i]])
		}

	if(single.cds && length(tmp.to)>1){ #combine multiple CDS in a cunning way
		#assume Exonerate produces the best cds first
		best.from <- tmp.from[[1]]
		best.to <- tmp.to[[1]]
		for(cdsi in 2:length(tmp.to)){ # run through remaining cds
			
			#test whether start/end of nth cds is within best.cds
			cds.overlap <- sapply(tmp.from[[cdsi]],function(X,best.from,best.to){any(X>best.from&X<best.to)},best.from,best.to)
			cds.overlap <- cds.overlap | sapply(tmp.to[[cdsi]],function(X,best.from,best.to){any(X>best.from&X<best.to)},best.from,best.to)
			
			#need to add lines to test whether start/end of best.cds is within nth cds
			#tmp.from.cdsi <- tmp.from[[cdsi]] #don't know why I need these... but crashes otherwise
			#tmp.to.cdsi <- tmp.to[[cdsi]]
			tmp.overlap <- sapply(tmp.from[[cdsi]],function(X,best.from){any(X<=best.from)},best.from) &
				sapply(tmp.to[[cdsi]],function(X,best.from){any(X>best.from)},best.from)
			cds.overlap <- cds.overlap | tmp.overlap
			tmp.overlap <- sapply(tmp.to[[cdsi]],function(X,best.to){any(X>=best.to)},best.to) &
				sapply(tmp.from[[cdsi]],function(X,best.to){any(X<best.to)},best.to)
			cds.overlap <- cds.overlap | tmp.overlap
			
			
			
			  #write separate sapply functions
			#cds.overlap <- cds.overlap | sapply(tmp.to[[cdsi]],function(X,best.to,tmp.from.cdsi){any(X>best.to&tmp.from.cdsi<best.to)},best.from,tmp.to.cdsi)
			#sapply(best.from,function(X,tmp.from.cdsi,tmp.to.cdsi){any(X>tmp.from.cdsi&X<tmp.to.cdsi)},tmp.from.cdsi,tmp.to.cdsi)
			#cds.overlap <- cds.overlap | sapply(best.to,function(X,tmp.from.cdsi,tmp.to.cdsi){any(X>tmp.from.cdsi&X<tmp.to.cdsi)},tmp.from.cdsi,tmp.to.cdsi)

			
			best.from <- c(best.from, tmp.from[[cdsi]][!cds.overlap])
			best.to <- c(best.to, tmp.to[[cdsi]][!cds.overlap])
			}
		tmp.from <- list(sort(best.from))
		tmp.to <- list(sort(best.to))
		}
	
	cds.old <- character(length=length(tmp.to))
	for(cds in 1:length(tmp.to)){
		for(j in 1:length(tmp.to[[cds]])){
			intron <- substr(ref.seq,start=tmp.from[[cds]][j],stop=tmp.to[[cds]][j])
			
			#get rid of codon cotaining N's from coding sequence
			#while(grep('^N',intron)) intron <- substr(intron,4,nchar(intron))
			#while(grep('N$',intron)) intron <- substr(intron,1,nchar(intron)-3)
			cds.old[cds] <- paste(cds.old[cds],intron,sep="")
			}
		
		} 
		
	#trim out masked regions
	cds.old <- gsub("#","",cds.old)

	#trim out codons with N's
	for(i in 1:length(tmp.to)){
		tmp <- character(1)
		for(codon.j in seq(1,nchar(cds.old[i]),3)){
			codon.nucl <- substr(cds.old[i],codon.j,codon.j+2)
			if(length(grep('N|n',codon.nucl))==0) tmp <- paste(tmp,codon.nucl,sep="")
			}
		cds.old[i] <- tmp
		}
	if(complementCDS){
		cds.old <- as.character(reverseComplement(DNAStringSet(cds.old)))
		}
	
	
	#translate, but use GeneR function to allow Ns
	old.prot <- as.character(sapply(cds.old,GeneR::strTranslate))



	list(transcript=cds.old,protein=old.prot)
	}

debug(getCDS)
tst.cds <- getCDS(seq.name=grouse600annot[2,"grouse.id"])
tst.cds <- getCDS(seq.name="cDNA_1788-1")


writeLines(muscle.char,"/Paterson/Datafiles/grouse/test_files/cDNA_1788-1_muscle.fa")

ta.char <- character(4)
ta.char[1] <- paste(">",grouse600annot$ensembl.id[grouse600annot$grouse.id=="cDNA_1788-1"],sep="")
ta.char[2] <- as.character(chick.seqs[[grep(tmp.chick.trans,names(chick.seqs))]])

ta.char[3] <- paste(">","cDNA_1788-1",sep="")
ta.char[4] <- tst.cds$transcript[1]
writeLines(ta.char,"/Paterson/Datafiles/grouse/test_files/cDNA_1788-1_ta.fa")

library(ape)

prot.afa <- read.BStringSet("/Paterson/Datafiles/grouse/test_files/cDNA_1788-1_muscle.afa")
dna.afa <- read.BStringSet("/Paterson/Datafiles/grouse/test_files/cDNA_1788-1_ta.afa")


fixedSubs <- function(chick.transcript,grouse.transcript,dir="/Paterson/Datafiles/grouse/test_files",transition.ratio=FALSE){
	#calculates fixed differences between a pair of sequences
	#(called chicken and grouse)
	require(GeneR)
	muscle.char <- character(4)
	muscle.char[1] <- paste(">ChickenProtein")
	muscle.char[2] <- strTranslate(chick.transcript)
	muscle.char[3] <- paste(">GrouseProtein")
	muscle.char[4] <- strTranslate(grouse.transcript)
	
	muscle.file <- paste(dir,"/proteins_tmpfile.fa",sep="")
	muscle.file.out <- paste(dir,"/proteins_tmpfile.afa",sep="")
	writeLines(muscle.char,muscle.file)
	
	ta.char <- character(4)
	ta.char[1] <- paste(">ChickenTranscript")
	ta.char[2] <- chick.transcript
	ta.char[3] <- paste(">GrouseTranscript")
	ta.char[4] <- grouse.transcript
	
	tranalign.file <- paste(dir,"/transcripts_tmpfile.fa",sep="")
	tranalign.file.out <- paste(dir,"/transcripts_tmpfile.afa",sep="")
	writeLines(ta.char,tranalign.file)

	#need to set R up use same PATH as X11
	system(paste("/usr/local/bioinf/muscle/muscle -in",muscle.file,"-out",muscle.file.out))
	system(paste("/opt/local/bin/tranalign",tranalign.file,muscle.file.out,tranalign.file.out))
	require(Biostrings)
	prot.afa <- read.BStringSet(muscle.file.out)
	dna.afa <- read.BStringSet(tranalign.file.out)
	
	tst <- strsplit(as.character(prot.afa),"")
	rm.gaps <- unique(c(grep("-",tst[[1]]),grep("-",tst[[2]])))
	tst2 <- list(chick.prot = toupper(tst[[1]][-rm.gaps]), grouse.prot = toupper(tst[[2]][-rm.gaps]))
	
	#insert homology criteria here
	# .. nothing yet ..
	
	#yn00 estimates
	paml.estimates <- try(pamlCalc(dna.afa,dir))
	#note, if paml fails may still take results from existing files
	
	prot.diff <- sum(tst2$chick.prot!=tst2$grouse.prot)
	tst <- strsplit(as.character(dna.afa),"")
	rm.gaps <- unique(c(grep("-",tst[[1]]),grep("-",tst[[2]])))
	tst2 <- list(chick.dna = toupper(tst[[1]][-rm.gaps]), grouse.dna = toupper(tst[[2]][-rm.gaps]))
	if(any(!tst2$chick.dna %in% c("A","G","C","T"))) warning("unexpected base in chicken sequence\n")
	if(any(!tst2$grouse.dna %in% c("A","G","C","T"))) warning("unexpected base in grouse sequence\n")
	dna.diff <- sum(tst2$chick.dna!=tst2$grouse.dna)
	if(transition.ratio==TRUE){
		#count up number of transitions and tranversions
		tst2$chick.type <- sapply(tst2$chick.dna,function(X){
			if(X %in% c("A","G")) return("R") #purine
			if(X %in% c("C","T")) return("Y") #pyrimidine
			warning('unknown base in chicken sequence')
			return("N") #only if function gets this far
			})
		tst2$grouse.type <- sapply(tst2$grouse.dna,function(X){
			if(X %in% c("A","G")) return("R") #purine
			if(X %in% c("C","T")) return("Y") #pyrimidine
			warning('unknown base in grouse sequence')
			return("N") #only if function gets this far
			})
		no.transversions <- sum(tst2$chick.type!=tst2$grouse.type)
		no.transitions <- dna.diff - no.transversions
		tmp <- c(prot.diff,dna.diff-prot.diff,no.transitions,no.transversions)
		names(tmp) <- c("Dn","Ds","transitions","transversions")
		}else{
		tmp <- c(prot.diff,dna.diff-prot.diff)
		names(tmp) <- c("Dn","Ds")
		}
	if(class(paml.estimates)=="try-error"){
		tmp <- c(tmp,rep(-1,9))
		}else{
		tmp <- c(tmp,unlist(paml.estimates))
		}
	tmp
	}

#debug(fixedSubs)
#fixedSubs(as.character(chick.seqs[[grep(tmp.chick.trans,names(chick.seqs))]]),tst.cds$transcript[1])
#work out how many coding and non-coding changes
pamlCalc <- function(sequences,dir){
	require(Biostrings)
	#sequences is a DNAstringset object
	#write sequences into phylip sequential format
	paml.file <- paste(dir,"/paml_input.phy",sep="")
	if(length(sequences[[1]])!=length(sequences[[2]])) stop('sequences of unequal length\n')
	paml.char <- character(length=1+2*length(sequences))
	paml.char[1] <- paste(length(sequences),length(sequences[[1]]))
	for(i in 1:length(sequences)){
		paml.char[2*i] <- names(sequences)[i]
		paml.char[1+2*i] <- toupper(as.character(sequences[[i]]))
		}
	writeLines(paml.char,paml.file)
	#assumes that there are control files yn00.ctl and codeml.ctl in dir
	#outputs to paml_output.out
	dir.now <- getwd()
	setwd(dir)
	system("yn00") #run Yang and Neilsen 200
	paml.grep <- system("grep -A 2 \'kappa\' paml_output.out",intern=TRUE)[c(1,3)]
	#1     2       3       4        5     6       7        8     9    10    11
	#seq. seq.     S       N        t   kappa   omega     dN +- SE    dS +- SE
	paml.out <- as.numeric(strsplit(paml.grep[2]," ")[[1]][grep('^[0-9]',strsplit(paml.grep[2]," ")[[1]])])
	setwd(dir.now)
	list(dN=paml.out[8],dN.se=paml.out[9],dS=paml.out[10],dS.se=paml.out[11],
		kappa=paml.out[6],omega=paml.out[7],Ls=paml.out[3],Ln=paml.out[4],t=paml.out[5])
	
	
	}

